####################################################################################
###
###   Immigration, Public Housing and Support for the French National Front
###   Gloria Gennaro
###
###   Paper Figure 3
###
####################################################################################

library(dplyr)
library(ggplot2)
library(AER)
library(multiwayvcov)

################################################################################
# Set Up
################################################################################

# Load Working Directories
source(paste0(wd_main, '/2_code/00_working_directories.R'))

# Load data and clean
source(paste0(wd_code, '/02_data_load_clean_rdd.R'))

# # Keep only years where there is an election
df = df_did[which(!is.na(df_did$fn)),]

# Keep only year>=2015, where I have demand data
df = df[!is.na(df$hlm_active_dem),]


################################################################################
# Variable definitions
################################################################################

df$imm_quant = cut(df$imm_share_99,
                   breaks=c(quantile(df$imm_share_99, probs = seq(0, 1, by = 1/3))),
                   labels=c("1","2","3"), include.lowest=T)

df$hlm_fr_rate = (1 + df$hlm_active_dem_fr)/(df$hlm_active_dem_fr + df$hlm_closed_dem_fr + 1)
df$hlm_exc_rate = (1 + df$hlm_active_dem_exc)/(df$hlm_active_dem_exc + df$hlm_closed_dem_exc + 1)


################################################################################
# Functions
################################################################################


LinearCombTest = function (lmObject, vars, .vcov = NULL) {
  if (is.null(.vcov)) .vcov = vcov(lmObject)
  beta = coef(lmObject)
  sumvars = sum(beta[vars])
  se = sum(.vcov[vars, vars]) ^ 0.5
  tscore = sumvars / se
  pvalue = 2 * pt(abs(tscore), lmObject$df.residual, lower.tail = FALSE)
  matrix(c(sumvars, se, tscore, pvalue), nrow = 1L,
         dimnames = list(paste0(paste0(vars, collapse = " + "), " = 0"),
                         c("Estimate", "Std. Error", "t value", "Pr(>|t|)")))
}


estimate_by_imm = function(outvar, bdw){

  temp = df[which(abs(df$running)<as.numeric(bdw)),]
  container = data.frame()
  
  print(paste0('############ ', outvar, ' ############'))

  controls = " + scale(hlm_share_99) + scale(res_owned_share_99) + scale(imm_share_99) + scale(fn_1995)"
  formula2sls = paste0(outvar, ' ~ sruT*factor(imm_quant) + running_pre', controls, ' | dummy*factor(imm_quant) + running_pre', controls)

  print(paste0('----- 2sls ----------'))
  model = ivreg(data=temp, as.formula(formula2sls))
  vcv = cluster.vcov(model, cluster = temp$agglo_name)
  model2 = coeftest(model, vcov = vcv)
  
  container = rbind(container, c(LinearCombTest(model, c("sruT"), vcv), 'Low Immigration', bdw ))
  container = rbind(container, c(LinearCombTest(model, c("sruT","sruT:factor(imm_quant)2"), vcv), 'Medium Immigration', bdw ))
  container = rbind(container, c(LinearCombTest(model, c("sruT","sruT:factor(imm_quant)3"), vcv), 'High Immigration', bdw ))
  
  colnames(container) = c('coef', 'se', 'tval', 'pval', 'group', 'bdw')
  container[,1:4] =  sapply(container[,1:4], function(x) round(as.numeric(as.character(x)), 3))
  
  output = list()
  output[[1]] = model
  output[[2]] = model2
  output[[3]] = container
  
  return(output)
}



##################################################################################
# Estimations
##################################################################################

contex1 = estimate_by_imm('scale(hlm_exc_rate)', bdw = 900)
contfr1 = estimate_by_imm('scale(hlm_fr_rate)', bdw = 900)


##############################################################################
#  Plot 
##############################################################################

container = rbind(contfr1[[3]], contex1[[3]])
container$variable = c(rep('Native\nUnmet Demand Rate', 3), rep('Immigrant\nUnmet Demand Rate', 3))
container$group.f = factor(container$group, levels = rev(unique(container$group)))
container$variable.f = factor(container$variable, levels = unique(container$variable))
container$label = paste0(as.character(round(container$coef,3)), '\n (',  as.character(round(container$se,3)), ')')

plot1 = ggplot(container, aes(x=coef, y=group.f)) + 
  geom_point(aes(x=coef, y=group.f), size=2) + 
  geom_errorbar(aes(xmax = coef+se*qnorm(.05,lower.tail=FALSE), xmin = coef-se*qnorm(.05,lower.tail=FALSE)), width =.2, size=.5)  +
  theme_bw() +
  xlab('') + 
  ylab('') +
  ggtitle("") +
  geom_text(aes(x = coef, y = group.f, label =label), size=4, vjust = 0, nudge_y = 0.1) +
  geom_vline(xintercept = 0, colour = gray(1/2), lty = 2) +
  facet_wrap(~variable.f, scales = "free_x") +
  theme(text = element_text(size=16),
        panel.grid.major.x = element_blank(),
        axis.text.y = element_text(size=12, hjust = 0.5))


plot1
ggsave(paste0(wd_res, '/figures/fig3.pdf'), height=7, width=10)





